Semiclassical approximation solved by Monte Carlo as an efficient impurity solver for 
dynamical mean field theory and its cluster extensions 
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We propose that a combination of the semiclassical approximation with Monte Carlo simulations 
can be an efficient and reliable impurity solver for dynamical mean field theory equations and their 
cluster extensions with large cluster sizes. In order to show the reliability of the method, we consider 
two test cases: (i) the single-band Hubbard model within the dynamical cluster approximation with 
4- and 8-site clusters and (ii) the anisotropic two-orbital Hubbard model with orbitals of different 
band width within the single-site dynamical mean field theory. We compare our results with those 
obtained from solving the dynamical mean field equations with continuous time and determinant 
quantum Monte Carlo. In both test cases we observe reasonable values of the metal-insulator critical 
interaction strength U c /t and, while some details of the spectral functions cannot be captured by 
the semiclassical approximation due to the freezing of dynamical fluctuations, the main features are 
reproduced by the approach. 

PACS numbers: 71.10.Fd,71.30.+h,71.27.+a 



I. INTRODUCTION 

The single-site dynamical mean field theory (DMFT) 
approach has been extensively employed to explore the 
properties of the Hubbard model and, in general, of 
strongly correlated materials 1 ^ 3 . Though the metal- Mott 
insulator transition can be successfully accounted for 
within the DMFT approximation where dynamical fluc- 
tuations alone are emphasized, interesting physical phe- 
nomena such as spin density waves and superconductiv- 
ity can not be properly described due to the absence of 
spatial fluctuations. Cluster-extensions of DMFT like 
cellular-dynamical mean field theory (CDMFT) 3 and dy- 
namical cluster approximation (DCA)« can take into 
account the inter-site spatial fluctuations within the size 
of the cluster in addition to the dynamical fluctuations. 
For example, cluster extensions of DMFT with 4-site 
clusters in combination with continuous-time quantum 
Monte Carlo (CT-QMC) or exact diagonalization (ED) 
can partly capture the physics of Fermi liquid (FL), non- 
FL (or pseudo gap), Mott insulator, and superconduc- 
tivity on an equal footing^—. However, due to the com- 
putational expense of CT-QMC and ED, the system size 
is still limited to only small clusters and not all physical 
properties can be equally precisely studied. In particular, 
the hybridization expansion CT-QMC approach is able 
to treat small cluster sizes up to N c = 4 due to an expo- 
nential increase of the local Hilbert space with TV" , . 13 ' 14 . 
The ED approach encounters a similar problem as the hy- 
bridization expansion CT-QMC does. Even though the 
determinant QMC (DQMC) and interaction expansion 
CT-QMC approaches are applicable to large-cluster sys- 
tems, the computational expense is proportional to the 
square of the interaction strength U, the temperature T, 



and the number of cluster sites N<M-~—. Therefore, a 
fast and reliable impurity solver for DMFT calculations 
and its cluster extensions with large cluster sizes is highly 
desirable. 

The semiclassical approximation (SCA) has been pro- 
posed as an impurity solver for DMFT and its cluster 
extensions^—. It gives reasonable results in the high 
temperature and strong interaction regions. The Fermi 
liquid behavior in the weak- coupling region is replaced 
by an anomalous metallic state due to the freezing of 
quantum fluctuations. Nevertheless, a big advantage of 
this solver is that self-energy and spectral functions can 
be directly obtained on the real frequency axis. This 
avoids the uncertainty from analytic continuation which 
has to be done in various QMC approaches. The pre- 
viously used SCA approach^— was limited though to 
consider small cluster sizes of N c = 4 due to the diffi- 
culty of the multi-dimensional integrations. 

In the present work, we propose to combine the SCA 
approach with the Monte Carlo (MC) method. The lat- 
ter is used to evaluate the multi-dimensional integrals. 
In this case the computational expense of the SCA ap- 
proach is only dependent on the cluster size N c and Mat- 
subara frequency, unlike the QMC approaches with de- 
pendence on U, T, N c , and Matsubara frequency. We 
apply our scheme to two test cases: (i) the one-orbital 
Hubbard model on the square lattice at half-filling within 
the DCA with cluster sizes N c = 4 and N c = 8 and (ii) 
the anisotropic two-orbital Hubbard model with differ- 
ent band widths on the Bethe lattice at half-filling within 
single-site DMFT. We benchmark our results with those 
obtained from solving the DCA and DMFT equations 
with CT-QMC and DQMC, respectively. For case (i) 
we find that even though the Fermi liquid behavior is 
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not obtained^, the critical onsite Coulomb interaction 
U c /t for the metal- insulator transition calculated by our 
SCA approach for both cluster sizes shows a reasonable 
agreement with the value obtained from CT-QMC which 
should be numerically exact. In particular, we find that 
the behavior of the density of states at the Fermi level in 
each momentum sector obtained from both SCA and CT- 
QMC approaches are quantitatively consistent with each 
other. For case (ii) we also find a reasonable agreement 
of U c /t obtained from SCA and DQMC. The orbital- 
selective phase transition is also correctly detected. How- 
ever, we observe that if the band- width difference be- 
tween narrow and wide orbitals is large, a causality prob- 
lem appears in the SCA results. 

The paper is organized as follows. In Sec. [IT] we present 
the general formalism of the semiclassical approximation 
and its application to cases (i) and (ii). In Sec. lIIIl we dis- 
cuss our SCA calculations and compare them with both 
CT-QMC and DQMC results and finally, in Sec. [TV] we 
summarize our findings. 



located at a distance £-th apart. The orthogonality is 
imposed by 
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j) n ii( T ) m t ne decoupling scheme can be written as: 
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where N% = (n^ + n,j.) and Mj = [n^ — n^) are the par- 
ticle number and magnetization, respectively. In terms 
of these definitions and within the SCA approximation, 
the partition function transforms into 

D [c+c] &tf dT ^ <^'V>i(T)MT,r')<Mr') + f SS M ?( T \ 

(6) 

where the Nf(r) term, which describes charge fluctua- 
tions, is neglected in the SCA approximation. This ex- 
pression can be rewritten as 



II. SEMICLASSICAL APPROXIMATION 



z= i d y c 



(7) 



In this section we will review the formalism of the SCA 
approach^— adapted to the two test cases considered 
in this work; the 8-site DCA and two-orbital DMFT sys- 
tems with paramagnetic solutions. 



A. General formalism 

The partition function can be written as: 

Z = J D[Jc]e- {So+s ^\ 



where 
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and 



dr / dr'^Ur)&a(r,T')Mr'), (3) 



So = 



where ip^ = (c\ a . . . c^) and c| CT (ci a ) is a Grassmann 
number corresponding to the Fermionic creation (annihi- 
lation) operator at site i and spin er, and a a = a (,aK^ a 
where a^ a are inverted frequency- dependent Weiss fields 
and K^c are n x n matrices defined according to the 
chosen cluster. Here, n denotes the number of sites in 
the multi-site system (case (i)) or two times the number 
of orbitals in the multi-orbital system (case(ii)), and £ 
denotes the distance between two sites within the clus- 
ter. For example, a^ 1 means the local (on-site) Weiss 
field while aT 1 is the inter-site Weiss field with the sites 



r^(r)a CT (^0^(^)-(|-^)). 



with 

rP 

S= dr 
Jo 

Here we assume that the new auxiliary fields 4>i (r) , which 
are given by a continuous Hubbard-Stratonovich trans- 
formation, arc t- independent {4>i{T) = (f>i). 
We replace M^(t) by 



M, 



(r)= f dr'J2cUr)af5(r-r')c ls ,(r'\ (9) 
Jo 



where a z is the third Pauli matrix. Via a Grassmann in- 
tegration and Fourier transformation, the partition func- 
tion (Eq. ) is finally given as 



AU 



hiDet 



-/3(a (T (iw„) + i^iffz) 



(10) 

where to n are Fermionic Matsubara frequencies and i is 
the dimension of integration. 

The impurity Green's function can be obtained as 
G^P^Wn) = G irJ (iw n )K ic r where G irT (iuj n ) is 



1 dlnZ 



n da^ a {ioj r , 



(11) 



where n is the normalization factor which is given by 
Eq. ((4]). The Green's function on the real frequency w is 
also calculated by Eq. (fTTj) with replacement of w„ into 
u+iS. In our calculations we consider a broadening factor 
<5 = 0.003. The i-dimension integration for classical fields 
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4>i is evaluated by the MC approach, where the weight 
function W (</>,) for MC simulations is given as 



(T,UJ n 
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(12) 



B. 8-site dynamical cluster approximation 

The partition function in the SCA approach is de- 
scribed in a real-space basis in Eqs. (JT|)-(l3j). The 8x8 
matrices of inversed Weiss fields (Eq. (fTU)) and Eq. (JTTJ)) 
in the 8-site DCA calculations for the Hubbard model on 
the square lattice are given as 
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where spin indices are omitted for simplicity and the nor- 
malization factors \ and are introduced in order to 

fulfill the orthogonality condition Eq. (((!])). The indices 
£ = 0,1,2, and 3 indicate on-site, 1st neighbor, 2nd 
neighbor, and 3rd neighbor, respectively. The cluster 
we used for constructing the a(iuj n ) matrices with pe- 
riodic boundary condition is shown in the left panel of 
Fig. [1] and the division of BZ for DCA calculations is 
presented in the right panel of Fig. [TJ The real-space 
impurity Green's functions in Eq. ( (jilt ) are more clearly 
expressed as 
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Det A(a 5 (iw n )) 
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where A(a^(iuj n )) = —/3(a(iuj n ) + \4>i<J z ) are 8x8 ma- 
trices. The impurity Green's function in Eq. ([Ft)) is 
measured by 
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(a) real space 






O 


8 
O 




O 


4 




r 


O 








3 






5 8 
O O 


/2 4 
/ • • 


5\ 
\ 


8 
O 


2 
O 


7 1 \ 
O O 


3 


6 

O 


7 


/ 
O 


3 
O 




5\ 


8 


/ 






3 6 
O O 


7 
O 


O 


3 
O 


6 

O 


7 
O 


(b) k space , 
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FIG. 1: (Color online) Cartoons for (a) the eight-site dynami- 
cal cluster approximation in real space and for (b) the division 
of Brilliouin zone in momentum space we used. 



where ij £ {1, 2} denote orbital indexes. U and U' are, re- 
spectively, onsite intra-orbital and inter-orbital Coulomb 
interaction parameters and J z is the Ising Hund's cou- 
pling term. We are not considering the spin-flip and pair 
hopping terms in our calculations. The inversed Weiss 
field is given as 
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We now decouple the interaction term Eq. (I16p using 
Eq. ©: 



C. 2-orbital dynamical mean field theory 



The interaction term 
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Neglecting charge fluctuations N? , the partition function 
can be written as 



fP TJ J J' 

Z__ j (j(Ml + Mi) + -(M : 



where U" = U'-J, and 
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where ip^ = (c[^,c\^,c\^,c\^) (compare with Eq. 
In order to make integration feasible, Eq. (|18[) is trans- 
formed into 
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FIG. 2: (Color online) The density of states at the Fermi 
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where we used the continuous Hubbard-Stratonovich 
transformation as in Eq. flSJ. Next (see Eq. ©), M^(t) 
is replaced by 



dr— 4>(,M^{t) level p(ui = 0) » /3G(§) and p(cj = 0), directly measured in 
* real-frequency space within DCA(SCA), as a function of U/t 

for N c = 4 and 8 at T/t = 1/12. We find the critical metal- 
insulator interactions U c /t = 4.1 and 5.2 for N c — 4 and 8, 
respectively. The error bars are smaller than symbol sizes. 
The deviations of between /3G(§) and p(u) = 0) are around 
ten percent. 
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Site indices are omitted due to the single-site DMFT cal- 
culation and ct| are 4x4 matrices: 
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Finally, the partition function is rewritten as 



where V(<j){) 
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(22) 

n is the number of 



4U' 1 AU 

orbitals (compare with Eq. (|10[) b The impurity Green's 
functions are calculated by Eq. (fT4|) . 



D. Monte Carlo measurement 

The weight functions for MC calculations have been 
given in Eq. (TT2l . We employ about 400 Matsubara fre- 
quencies for performing the frequency sum in Eq. (1121) . 



The number of classical holds 4>i is the same as the num- 
ber of cluster sites in the DC A calculations (case (i)). 
For case (ii) the number of classical fields is given by 
(2m — 1) x to, where to is the number of orbitals. In 
order to avoid a local minimum problem in the MC cal- 
culation, we use around forty initial configurations and 
we perform about 4 x 10 5 MC samplings for each different 
initial configuration. 



III. RESULTS 

A. The metal-insulator transition in the 8-site 
dynamical cluster approximation 

In what follows, we will show the reliability of our SCA 
impurity solver by presenting the results obtained for 4- 
and 8- site DCA calculations for a two-dimensional Hub- 
bard model on the square lattice at half- filling (case (i)). 

First, we compare the critical U c /t obtained from 
DCA(SCA) calculations with that obtained from 
DCA(CT-QMC) with N c = 4 and 8. Within DCA(CT- 
QMC)2&24 U c /t = 4.5 (N c = 4) and U c /t = 6.5 (N c = 8). 
We would like to note that for the two-dimensional Hub- 
bard model on the square lattice at half-filling, the Slater 
gap (band gap) should be opened for an infinitesimal 
U due to long-range antiferromagnetic spatial correla- 
tions^^. This suggests that we should expect a smaller 
U c /t for increasing cluster sizes. However, this is not 
what it is observed above. We think that the reason for 
this discrepancy lies on the fact that plaquette singlet 
ordered states become more favorable for N c = 4 than 
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for N c = 8 and artificially stabilize an insulating state 
in N c — 4. We now check whether this behavior is also 
observed in the SCA approach. In Fig. [2] we plot the 
DCA(SCA) density of states at the Fermi level obtained 
as (1) p(u = 0) fa /3G(f) and (2) directly calculated 
in real- frequency space p(u = 0), as a function of U/t 
for N c = 4 and 8 at T/t = 1/12. The imaginary time 
Green's function G(t) is calculated by the Fourier trans- 
formation of G(iw n ) in Eq. (fTTI) . We find U c /t — 4.1 and 
5.2 for N c — 4 and 8 respectively. The trend, i.e., smaller 
critical interaction U c /t for N c = 4 than for iV c = 8, is 
the same as in DCA(CT-QMC). We also detect that the 
critical interactions U c /t in the SCA method are smaller 
than those calculated by the CT-QMC method. The rea- 
son why the insulating state is overestimated, is due to 
the fact that the auxiliary field is assumed to be r in- 
dependent in SCA, indicating a freezing of dynamical 
fluctuations in SCA. 

Usually one always employs the relation p(uj = 0) ps 
/8Gr(%) to determine the critical interaction U c /t for the 
metal-insulator transition in CT-QMC. This is done in 
order to avoid the performance of an analytical continua- 
tion, which will introduce some uncertainties. Therefore, 
it is interesting to check whether this relation is valid in 
all cases. Since DCA(SCA) provides results directly in 
real- frequency space, both definitions can be tested on 
the same footing. Fig. [5] shows that both definitions 
are in good agreement in the weak-coupling and strong- 
coupling regions, while they show deviations of about ten 
percent close to the critical value. In addition, we find 
a causality problem in the weak-coupling region (for U/t 
values smaller than 3.0). 

Next, we would like to analyze the spectral functions in 
different DCA cluster momentum sectors, i.e., A{K.,u) 
at K = (0, 0), (0, tt), (f , | ) and (tt, tt) shown in Fig. CD (b) 
for several values of U/t at T/t = 1/12. In Figs. [3] (a) 
we show the non- interacting case (U/t — 0.0). While 
the spectral functions A(K,lu) at K = (0,0) and (%, tt) 
sectors are well separated from each other, resembling 
the behavior of band insulators, the spectral functions 
at K = (0, 7r) and (§,§) sectors cross the Fermi level, 
showing metallic behavior. The van-Hove singularity is 
present in the spectral function in the K = (tt, 0) / (0, tt) 
sector. The behavior of A(K, ui) for N c — 8 is comparable 
to N c — 4 results^ and can be understood in terms of the 
non-interacting band structure. In Figs. [3] (b) we show 
A(K, u) for U/t = 3.2 (weak-coupling region). A(K,w) 
at K = (0, 0) and (tt, tt) intersect with each other due to 
the band splitting and spectral weight transfer induced 
by U. This is an indication of Mott physics. The insu- 
lating behavior still remains in these two sectors and the 
whole band width is slightly narrowed due to correlation 
effects. The van-Hove singularity that was present in the 
K = (tt, 0) sector in the non-interacting case, is dramat- 
ically suppressed with increasing U /t. The absence of a 
strong quasi-particle peak in the weak interaction region 
is due to the freezing of dynamical fluctuations. This is 
a shortcoming of the SCA method. As the interaction 




FIG. 3: (Color online) The spectral functions A(K, cj) in he 
different DCA momentum sectors K for (a) U/t = 0.0, (b) 
3.2 (c) 4.0 and (d) 6.0 at T/t = 1/12. 



U/t is increased, a pseudo-gap behavior is present with 
suppression of the spectral functions at K = (0, tt) and 
(§, f ) at U/t = 4.0 (Fig. 1 (c)), and the Mott insula- 
tor appears in the strong-coupling region at U/t = 6.0 
(Fig. |3](d)). 

In Fig. [4] we show /3G(K, |) as a function of U/t at 
T/t = 1/12. When the interaction U/t is turned on, the 
values of /3G(K, f ) at K = (0, 0) increase until U/t = 3.6 
due to spectral weight transfer caused by electronic corre- 
lations, and a metallic behavior is seen in the K = (0, 0) 
sector in the intermediate interaction strength regions. 
When the interaction becomes strong, the gap opens 
and /3G(K, ^) goes to zero. In the K = (0, tt) sector 

/3G(K, f) decreases monotonously and goes to zero at 
U/t = 5.6. In the K = (■§,§) sector, the values of 
/?G(K, remain nearly constant up to U/t = 3.0. Be- 
yond U/t = 3.0, they decrease and the gap opens com- 
pletely around U/t = 5.6. 

Finally, we compare the results in Fig.|4]to those in Fig. 
8 of Ref. [13 calculated within the CT-QMC approach. 
The behavior at K = (0,tt) in both SCA and CT-QMC 
approaches is qualitatively the same, even though the 
critical interactions U c /t are different. The main differ- 
ence between these two results is in the K = ( § , f ) sec- 
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FIG. 4: (Color online) The quantity /3G(K, f ) in the different 
DCA momentum sectors K in the 8-site DCA approach at 
T/t = 1/12. The values of /3G(K, f ) at K = (tt, tt) are the 
same as those at K = (0,0). The error bars are smaller 
than symbol sizes. These results are comparable with those 
obtained from the interaction-expansion CT-QMC approach 
in Ref. [M 



tor. The CT-QMC results indicate a first-order transition 
with a discontinous behavior of /3G(K, at the critical 
interaction U c /t, while the SCA results show a continu- 
ous transition with a smooth decrease of j3G(K, ^). We 
think that this discrepancy between the two approaches 
also comes from the approximation that the dynamical 
fluctuations are frozen in the SCA method. 



B. Orbital-selective phase transitions in the two- 
orbital dynamical mean field theory 

The orbital-selective phase transition (OSPT), where 
metallic behavior is seen in the wide band while a metal- 
insulator transition is observed in narrow band, has been 
intensively studied in model systems as well as real ma- 
terials during the last ten years22r— . We study the 
anisotropic two-orbital Hubbard model with a narrow 
bandwidth of Wi = 2.56 (ti = 0.8) for the first orbital 
and a wide bandwidth of W2 — 4 (t 2 = 1.0) for the 
second orbital at half-filling on the Bethe lattice using 
single-site DMFT. The interaction part of the Hamilto- 
nian is given by Eq. with U' = f and J z = j. The 
DOS results in both orbitals iobtained with the SCA are 
shown in Fig.[5](a)-(c). Metallic behaviour in both bands 
is observed for the weak coupling strength U/t2 = 2.0 
in Fig. [S] (a). As the interaction U/t is increased, the 
OSPT behavior is present in the intermediate regions for 
U/t 2 = 2.8 (see Fig. |S](b)). Finally, insulating states in 
both orbitals are seen in the strong-coupling region for 
U/t 2 = 3.6 in Fig. m (c). 

In the following we compare the critical values of the 
interaction strength obtained from the SCA and DQMC 



calculations in the case of bandwidths of W\ = 2 and 
W2 = 4 in the narrow and wide bands, respectively. 
The critical values in the DQMC approach are given as 
U Cl /t2 = 2.0 and U C2 /t2 = 2.8 in narrow and wide bands, 
respectively^. Our SCA results show the critical val- 
ues U Cl /t2 = 2.0 and U C2 /t2 = 3.2 from the analysis of 
the Green's function in the Matsubara frequency space 
in good agreement with DQMC. We encounter though 
a causality problem if the difference between the band 
widths of narrow and wide orbitals is large. 



IV. SUMMARY 

In summary, in this paper we propose that the semi- 
classical approximation in combination with the Monte 
Carlo method can be used to study large clusters and 
multi-orbital systems and is easy to embed into DMFT 
and its cluster extensions. We investigate the single- 
orbital Hubbard model by the DCA(SCA) method with 
cluster sizes of iV c = 4 and 8 and a two- orbital sys- 
tem by DMFT(SCA). The results are compared with 
existing DCA(CT-QMC) and DMFT(CT-QMC) results. 
The critical interactions U/t of SCA and CT-QMC ap- 
proaches in both cases are in reasonable agreement. In 
the 8-site DCA cluster calculation, we analyze the spec- 
tral functions A(K,cj) at each momentum sector. The 
only difference between SCA and CT-QMC results is that 
the CT-QMC shows a discontinuous behavior of the spec- 
tral density in the K = (^, ^) sector around the critical 
value of interactions, while the SCA exhibits smoothly 
decreasing behavior. We think that the reason for the 
discrepancy is that quantum fluctuations are frozen in 
the SCA approach. In the two-orbital DMFT(SCA) cal- 
culation, we observe the orbital selective phase transition 
as in previous studies. This method is rather powerful 
since it can be applied for problems where other impu- 
rity solvers remain computationally unfeasable but one 
should be aware of possible causality problems in some 
cases. 
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FIG. 5: (Color online) The density of states p(oj) for (a) U/t 2 = 2.0, (b) 2.8 and (c) 3.6 at T/t 2 = 1/12.0. The bandwidths 
for narrow and wide orbitals on the Bethe lattice are W\ = 2.56 (ti = 0.8) and W2 = 4.0 (£2 = 1-0), respectively. The metal, 
orbital-selective phase, and Mott insulator are present in (a), (b), and (c), respectively. 
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